#ifndef SPDE_ONE_MODE_H
#define SPDE_ONE_MODE_H

#include "Signal.h"

/**************************************************************
class SPDE_One_Mode solve the following 2D EW equation
\[
\frac{\partial h}{\partial t} = w(x,y)+c_2\left(\frac{\partial^2 h}{\partial x^2}+\frac{\partial^2 h}{\partial y^2}\right)+\xi(x,y,t)
\]
where $x\in[0,XYmax],y\in[0,XYmax]$. $\xi(x,y,t)$ is a Gaussian
white noise with zero mean and the following covariance
\[
\left<\xi(x,y,t)\xi(x',y',t')\right>=\sigma^2\delta(x-x')\delta(y-y')\delta(t-t').
\]
$w$ is a patterned deposition rate profile with the form
\[
w(x,y) = W+A\sin(\frac{2*k*\pi}{L}).
\]
Boundary condition:
h(0,y,t) = h(L,y,t),    h(x,0,t) = h(x,L,t)
\frac{\partial h}{\partial x}(0,y,t) = \frac{\partial h}{\partial x}(L,y,t)
\frac{\partial h}{\partial y}(x,0,t) = \frac{\partial h}{\partial y}(x,L,t)

Initial condition:
h(x,y,0) = h_0(x,y) = 0

In the code:
XYmax = L
Freq = k

Class SPDE_One_Mode use model decomposition method to
solvet the EW equation. Detailed can be found elsewhere.
**************************************************************/

class SPDE_One_Mode:public block{
  // Inputs =====================================
  double T;        // unit: K
  double W;        // unit: nm/s
  double A;        // unit: nm/s
  double Freq;     // unit: number of periods in the domain

  // Model parameters ===========================
  int    LatticeSize; // number of discrete points in the h profile
  double XYmax;    // unit: nm
  int    p;
  int    m;
  int    n;

  // Numerical solution parameters =============
  // For analytical expression
  double sqrt_dt;

  // States =====================================
  double *z;
  double *z_next;
  
  // Model parameters ===========================
  double lambda; // unit: 1/s

  // Random number generator ====================
  const gsl_rng_type * noise_type;
  gsl_rng * noise;

  // Fitted parameter model =====================
  double eq_c2(double T,double W);
  double eq_sigma2(double T,double W);

  // Eigenfunctions =============================
  inline double phi3mn(int m,int n,int ix,int iy){
    double coeff;
    if(n==0||n*2 == LatticeSize){
      coeff = SQRT_TWO/XYmax;
    }else{
      coeff = 2.0/XYmax;
    }
    return coeff*sin(2*m*pi/LatticeSize*ix)*cos(2*n*pi/LatticeSize*iy);
  };
 public:
  SPDE_One_Mode(int p,int m,int n,int L,double idt,double iXYmax);
  ~SPDE_One_Mode();
  virtual void update(double sysTime);
  virtual void reset();
  double getR2(); // Return r^2, unit = nm^2
  void setPara(int para_id,int para_val);
};
#endif
